##### SCRIPT QUE GERA A TABELA 3 DO APENDICE #####

# 1. CARREGA PACOTES ------------------------------------------------------

library(tidyverse)
library(janitor)
library(naniar)
library(mice)
library(xlsx)

# 2. CARREGA DADOS -------------------------------------------------------

#' Nesta etapa, foi implementada a leitura dos dados combinados de Enade-RAIS

# 3. IMPLEMENTA MODELOS DE REGRESSAO LOGISTICA ---------------------------

# 3.1 Modelo separado por área ----

imp_med <- filter(dados_graduados_empregados, area_curso_agreg == "Medicina")
imp_soc <- filter(dados_graduados_empregados, area_curso_agreg == "Ciências Sociais Aplicadas")
imp_ctem <- filter(dados_graduados_empregados, area_curso_agreg == "CTEM")

imp_dir <- filter(dados_graduados_empregados, area_curso_agreg == "Direito")
imp_edu <- filter(dados_graduados_empregados, area_curso_agreg == "Educação")
imp_eng <- filter(dados_graduados_empregados, area_curso_agreg == "Engenharia")

imp_hum <- filter(dados_graduados_empregados, area_curso_agreg == "Humanidades")
imp_sau <- filter(dados_graduados_empregados, area_curso_agreg == "Saúde e Bem-Estar")
imp_tec <- filter(dados_graduados_empregados, area_curso_agreg == "Tecnólogos")

f <- "ind_empregado ~ tp_sexo + raca_cor + idade_faixa + edu_fam + nt_ce_quartil_area + nt_fg_quartil + trab_grad + igc_faixa_4c + setor_ies_bin + regiao_curso"

model_med <- with(imp_med, glm(as.formula(f), family = binomial))
model_soc <- with(imp_soc, glm(as.formula(f), family = binomial))
model_ctem <- with(imp_ctem, glm(as.formula(f), family = binomial))
model_dir <- with(imp_dir, glm(as.formula(f), family = binomial))
model_edu <- with(imp_edu, glm(as.formula(f), family = binomial))
model_eng <- with(imp_eng, glm(as.formula(f), family = binomial))
model_hum <- with(imp_hum, glm(as.formula(f), family = binomial))
model_sau <- with(imp_sau, glm(as.formula(f), family = binomial))
model_tec <- with(imp_tec, glm(as.formula(f), family = binomial))

combine_med <- pool(model_med)
combine_soc <- pool(model_soc)
combine_ctem <- pool(model_ctem)
combine_dir <- pool(model_dir)
combine_edu <- pool(model_edu)
combine_eng <- pool(model_eng)
combine_hum <- pool(model_hum)
combine_sau <- pool(model_sau)
combine_tec <- pool(model_tec)

# 4. CALCULA PROBABILIDADES PREDITAS -------------------------------------

# 4.1 Por CE, area e setor da IES -----

# Funcao que cria dados-base para calculo de probabilidades preditas

prepara_new_data <- function(x = data, y) {
      x %>% 
            dplyr::distinct(area_curso_agreg) %>%
            filter(area_curso_agreg == y) %>% 
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            dplyr::arrange(area_curso_agreg) %>% 
            dplyr::mutate(tp_sexo = "Feminino",
                          raca_cor = "Negro/Indígena",
                          edu_fam = "Ensino Médio",
                          nt_ce_quartil_area = rep(c("q1", "q4"), 2),
                          nt_fg_quartil = "q3",
                          trab_grad = c("Não trabalha"),
                          igc_faixa_4c = "3",
                          idade_faixa = "25-",
                          setor_ies_bin = rep(c("Pública", "Privada"), each = 2),
                          regiao_curso = "SE") %>% 
            mutate_if(is.character, as.factor)
}

new_data_soc <- prepara_new_data(dados_graduados_empregados, "Ciências Sociais Aplicadas")
new_data_cte <- prepara_new_data(dados_graduados_empregados, "CTEM")
new_data_dir <- prepara_new_data(dados_graduados_empregados, "Direito")
new_data_edu <- prepara_new_data(dados_graduados_empregados, "Educação")
new_data_eng <- prepara_new_data(dados_graduados_empregados, "Engenharia")
new_data_hum <- prepara_new_data(dados_graduados_empregados, "Humanidades")
new_data_med <- prepara_new_data(dados_graduados_empregados, "Medicina")
new_data_sau <- prepara_new_data(dados_graduados_empregados, "Saúde e Bem-Estar")
new_data_tec <- prepara_new_data(dados_graduados_empregados, "Tecnólogos")

# Funcao que calcula probabilidades preditas

predict_prob <- function(x, y, z) {
      
      p <- x$analyses[[1]]
      p$coefficients = summary(y)$estimate
      predict(p, newdata = z,
              type = "response", se.fit = TRUE)
      
}

pred_soc <- predict_prob(model_soc, combine_soc, new_data_soc)
pred_cte <- predict_prob(model_ctem, combine_ctem, new_data_cte)
pred_dir <- predict_prob(model_dir, combine_dir, new_data_dir)
pred_edu <- predict_prob(model_edu, combine_edu, new_data_edu)
pred_eng <- predict_prob(model_eng, combine_eng, new_data_eng)
pred_hum <- predict_prob(model_hum, combine_hum, new_data_hum)
pred_med <- predict_prob(model_med, combine_med, new_data_med)
pred_sau <- predict_prob(model_sau, combine_sau, new_data_sau)
pred_tec <- predict_prob(model_tec, combine_tec, new_data_tec)

prepara_tabela <- function(x, y, z) {
      data.frame(z[,c("area_curso_agreg", "nt_ce_quartil_area", "setor_ies_bin")],
                 pred = x) %>%
            mutate(area_curso_agreg = y) %>% 
            select(1, 2, 3, 4) %>% 
            spread(nt_ce_quartil_area, pred.fit)
}

tab_pred_ce <- bind_rows(prepara_tabela(pred_soc, "SOC", new_data_soc),
                         prepara_tabela(pred_cte, "CTEM", new_data_cte),
                         prepara_tabela(pred_dir, "DIR", new_data_dir),
                         prepara_tabela(pred_edu, "EDU", new_data_edu),
                         prepara_tabela(pred_eng, "ENG", new_data_eng),
                         prepara_tabela(pred_hum, "HUM", new_data_hum),
                         prepara_tabela(pred_med, "MED", new_data_med),
                         prepara_tabela(pred_sau, "SAU", new_data_sau),
                         prepara_tabela(pred_tec, "TEC", new_data_tec)) %>% 
   data.frame

# 4.2. Por FG, area e setor ----

# Funcao que cria dados-base para calculo de probabilidades preditas

prepara_new_data <- function(x = data, y) {
      x %>% 
            dplyr::distinct(area_curso_agreg) %>%
            filter(area_curso_agreg == y) %>% 
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            dplyr::arrange(area_curso_agreg) %>% 
            dplyr::mutate(tp_sexo = "Feminino",
                          raca_cor = "Negro/Indígena",
                          edu_fam = "Ensino Médio",
                          nt_ce_quartil_area = "q3",
                          nt_fg_quartil = rep(c("q1", "q4"), 2),
                          igc_faixa_4c = "3",
                          trab_grad = c("Não trabalha"),
                          idade_faixa = "25-",
                          setor_ies_bin = rep(c("Pública", "Privada"), each = 2),
                          regiao_curso = "SE") %>% 
            mutate_if(is.character, as.factor)
}

new_data_soc <- prepara_new_data(dados_graduados_empregados, "Ciências Sociais Aplicadas")
new_data_cte <- prepara_new_data(dados_graduados_empregados, "CTEM")
new_data_dir <- prepara_new_data(dados_graduados_empregados, "Direito")
new_data_edu <- prepara_new_data(dados_graduados_empregados, "Educação")
new_data_eng <- prepara_new_data(dados_graduados_empregados, "Engenharia")
new_data_hum <- prepara_new_data(dados_graduados_empregados, "Humanidades")
new_data_med <- prepara_new_data(dados_graduados_empregados, "Medicina")
new_data_sau <- prepara_new_data(dados_graduados_empregados, "Saúde e Bem-Estar")
new_data_tec <- prepara_new_data(dados_graduados_empregados, "Tecnólogos")

# Funcao que calcula probabilidades preditas

predict_prob <- function(x, y, z) {
      
      p <- x$analyses[[1]]
      p$coefficients = summary(y)$estimate
      predict(p, newdata = z,
              type = "response", se.fit = TRUE)
      
}

pred_soc <- predict_prob(model_soc, combine_soc, new_data_soc)
pred_cte <- predict_prob(model_ctem, combine_ctem, new_data_cte)
pred_dir <- predict_prob(model_dir, combine_dir, new_data_dir)
pred_edu <- predict_prob(model_edu, combine_edu, new_data_edu)
pred_eng <- predict_prob(model_eng, combine_eng, new_data_eng)
pred_hum <- predict_prob(model_hum, combine_hum, new_data_hum)
pred_med <- predict_prob(model_med, combine_med, new_data_med)
pred_sau <- predict_prob(model_sau, combine_sau, new_data_sau)
pred_tec <- predict_prob(model_tec, combine_tec, new_data_tec)

prepara_tabela <- function(x, y, z) {
      data.frame(z[,c("area_curso_agreg", "nt_fg_quartil", "setor_ies_bin")],
                 pred = x) %>%
            mutate(area_curso_agreg = y) %>% 
            select(1, 2, 3, 4) %>% 
            spread(nt_fg_quartil, pred.fit)
}

tab_pred_fg <- bind_rows(prepara_tabela(pred_soc, "SOC", new_data_soc),
          prepara_tabela(pred_cte, "CTEM", new_data_cte),
          prepara_tabela(pred_dir, "DIR", new_data_dir),
          prepara_tabela(pred_edu, "EDU", new_data_edu),
          prepara_tabela(pred_eng, "ENG", new_data_eng),
          prepara_tabela(pred_hum, "HUM", new_data_hum),
          prepara_tabela(pred_med, "MED", new_data_med),
          prepara_tabela(pred_sau, "SAU", new_data_sau),
          prepara_tabela(pred_tec, "TEC", new_data_tec)) %>% 
   data.frame

# 5. EXPORTA TABELA -------------------------------------------------------

tab_pred_ce %>% 
   xlsx::write.xlsx(., "output/Tabela A3.xlsx", 
                    append = T, 
                    row.names = F, 
                    sheetName = "raw_ce")

tab_pred_fg %>% 
   xlsx::write.xlsx(., "output/Tabela A3.xlsx", 
                    append = T, 
                    row.names = F, 
                    sheetName = "raw_fg")
